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ABSTRACT 


We study variants of hierarchical modular network models suggested by Kaiser and Hilgetag [Front, in Neuroinform., 4 (2010) 
8] to model functional brain connectivity, using extensive simulations and quenched mean-field theory (QMF), focusing on 
structures with a connection probability that decays exponentially with the level index. Such networks can be embedded 
in two-dimensional Euclidean space. We explore the dynamic behavior of the contact process (CP) and threshold models 
on networks of this kind, including hierarchical trees. While in the small-world networks originally proposed to model brain 
connectivity, the topological heterogeneities are not strong enough to induce deviations from mean-field behavior, we show 
that a Griffiths phase can emerge under reduced connection probabilities, approaching the percolation threshold. In this 
case the topological dimension of the networks is finite, and extended regions of bursty, power-law dynamics are observed. 
Localization in the steady state is also shown via QMF. We investigate the effects of link asymmetry and coupling disorder, 
and show that localization can occur even in small-world networks with high connectivity in case of link disorder. 


Introduction 

In neuroscience, the criticality hypothesis asserts that the brain is in a critical state, at the boundary between sustained activity 
and an inactive regime. Theoretical and experimental studies show that critical systems exhibit optimal computational prop¬ 
erties, suggesting why criticality may have been selected in the evolution of the nervous system. * Although criticality has 
been observed in cell cultures,brain slices and anesthetized animals,the debate regarding criticality in alert animals and 
humans continues.®"^ Thus the criticality hypothesis remains controversial in brain science; for a review see.^ 

Normally, for a system to be at criticality, certain control parameters need to be tuned precisely, raising the perennial 
question of how such tuning is achieved in the absence of outside intervention. The possibility of self-tuning is well known in 
statistical physics; the paradigm of self-organized criticality (SOC) has been studied since the pioneering work of.'** Simple 
homogeneous models such as the stochastic sandpile exhibit criticality with power laws both in statics and dynamics. This 
has been understood as the result of a control mechanism that forces the system to an absorbing-state phase transition** 

Real nervous systems, however, are highly inhomogeneous, so that one must take into account the effects of heterogeneities. 
It is well known in statistical physics that disorder can cause rare-region (RR) effects*^ that smear the phase transitions. These 
effects can make a discontinuous transition continuous, or generate so-called Griffiths phases (GP),*^ in which critical-like 
power-law dynamics appears over an extended region around the critical point. In these regions, moreover, non-Markovian, 
bursty behavior can emerge as a consequence of a diverging correlation time. The inter-communication times of the nodes 
(which possess no internal memory) follow a fat-tailed distribution.*'* Thus, heterogeneities widen the critical region and can 
contribute to the occurrence of power laws. This provides an alternative mechanism for critical-like behavior without fine 
tuning, although attaining the GP does require some rough tuning, of the sort that is not difficult to find in biological systems. 

It was shown recently that the topological heterogeneity of the underlying networks can result in GPs in finite-dimensional 
systems*^ and can be a reason for the working memory in the brain.*® Although many networks exhibit the small-world 
property and so have an infinite topological dimension, naturally occurring networks are always finite, exhibit cutoffs, and 
therefore GPs can be expected as a consequence of inhomogeneous topology.*^ 

Many real networks can be partitioned seamlessly into a collection of modules. Each module is expected to perform an 
identifiable task, separate from the function of others.*^ It is believed that the brain is organized hierarchically into modular 
networks across many scales, starting from cellular circuits and cortical columns via nuclei or cortical areas to large-scale 
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units such as visual or sensory-motor cortex. At each level, connections within modules are denser than between different 
modules. Although empirical data confirm this modular organization on some scales,the detailed organization of brain 
networks is not yet experimentally accessible. 

Two particular kinds of hierarchical modular networks (HMN-l,HMN-2) model were proposed and investigated numer¬ 
ically and analytically in.^"* On large-world HMNs, which imply a finite topological dimension D, models of the spread of 
activity exhibit power-law dynamics and rare region effects. However, these power laws are system-size dependent, so that 
true GP behavior has not yet been proven. These authors also simulated spreading on empirical brain networks, such as the 
human connectome and the nervous system of C. elegans. Available empirical networks are much smaller than the synthetic 
ones and deviations from power laws are clearly visible. Both anatomical connections^^ and the synchronization networks of 
cortical neurons^^ suggest small-world topology.^® The brain network modules of^'* are weakly coupled in such a way that 
these HMNs are near the percolation threshold, as in the case of the models introduced in.'^ Note, however that requiring 
the network to be near percolation again raises tuning and robustness issues. Having weaker ties would lead to fragmented 
networks, while stronger ties result in infinite D and the absence of a GP. In the present work we do not assume such fine 
tuning; we maintain a high density of short edges, rendering the networks well connected. Another way of preserving the 
integrity of HMN networks with finite dimension involves the random tree-like structures studied here. 

To study synchronization,^^ commonly expected in brain dynamics, the Kuramoto model^^ has been implemented^*^ on 
the same networks as studied in.^"* In this case even weaker rare-region effects were found, resulting in stretched exponential 
dynamics. One of the main purposes of our study is to delineate conditions a GP. While the networks studied are of finite 
size, repeating the process on many network realizations and averaging over them, we clearly see convergence towards GP 
dynamics. We assume that multiple random network realizations may occur in the brain over time, due to reconfigurations 
of the synapses or as a consequence of weakly coupled sub-modules and changes in overall intensity of connections between 
different brain regions due to neuromodulation. 

In addition to the dimension, other topological features, which have yet to be studied systematically, may also influence 
RR effects. The clustering coefficient, defined as 

where n,- denotes the number of direct links interconnecting the k,- nearest neighbors of node /, exhibits different scaling 
behavior in modular networks than in unstructured networks. While in SF networks it decays as C{N) ~ and in 

random networks as C{N) ~ in HMNs it is constant. Furthermore, while in SF networks the clustering is degree-number 
independent, in HMNs it decays asC(k) withjS in the range 0.75 - 1.^' This means that the higher the node degree, the 

smaller is its clustering coefficient, possibly reducing its infectiousness in an epidemic process. This suggests an enhancement 
of localization as in SF models with disassortative weighting schemes. 

In this work we study spreading models defined on HMNs embedded in 2D space, and investigate how various inho¬ 
mogeneities (topological or intrinsic) contribute to the occurrence of localization and GPs. We begin by considering purely 
topological heterogeneity; later we study cases with disorder in the connections as well. It is well known that intrinsic link 
disorder can appear as the consequence of asymmetry, connecting nodes in the cortex (see^^). It has also been demonstrated 
that weights can vary over 4-6 orders of magnitude; most (including the majority of long-range connections) are actually quite 
weak.^'^ 

The relation of RR effects to localization in the steady state has been studied recently.*^ Here we discuss localization 
results obtained via a quenched mean-field (QMF) approximation. 

In addition to serving as models of brain connectivity, hierarchical modular networks are abundant in nature. They occur in 
social relations,in the WWW,^® metabolic^^ and information systems^^ for example. Thus the results reported here should 
find application in these and related areas. 

Hierarchical modular networks 

Recent studies indicate that activity in brain networks persists between the extremes of rapid extinction and global excitation. 
This so-called limited sustained activity (LSA) does not occur in spreading models defined on structureless, small-world 
networks. On the other hand, brain networks are known to exhibit hierarchical modular structure. This motivated Kaiser and 
Hilgetag (KH) to perform numerical studies on such networks, to investigate topological effects on LSA.^® Their hierarchical 
model reflects general features of brain connectivity on large and mesoscopic scales. Nodes in the model were intended to 
represent cortical columns rather than individual neurons. The connections between them were taken as excitatory, since there 
appear to be no long-distance inhibitory connections within the cerebral cortex.^® 

Kaiser and Hilgetag generated networks beginning at the highest level, and adding modules to the next lower level, with 
random connectivity within modules. They explored hierarchical networks with different numbers of levels, and numbers of 
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sub-modules at each level. The average degree (over all nodes in the network) was set to (k) = 50, motivated by experimental 
studies. They investigated different topologies by varying the edge density across the levels. All the networks studied by KH 
are of small-world type, i.e., they have an infinite topological dimension. 

The spreading model investigated by KH is a two-state threshold model, in which, at each time step, inactive nodes are 
activated with probability X if at least m of their neighbors are currently active, and active nodes deactivate spontaneously with 
probability v. This model is very similar to reaction-diffusion models known in statistical physics,^^^^ with a synchronous 
cellular automaton (SCA) updates. Starting from an initial state in which a localized set of nodes is active, these authors 
followed the density of active sites and their localization up to f ~ 200 time steps on networks of sizes A < 512. 

Kaiser and Hilgetag found that LSA can be found in a larger parameter range in HMNs, as compared with random and 
non-hierarchical small-world networks. The optimal range of LSA was found in networks in which the edge density increased 
from the top level {I = Imax) to the bottom (/ = li). Such topologies foster activity localization and rare-region effects. 

In this paper we investigate HMNs which possess increasing edge density from top to bottom levels, as did KH, but with 
finite topological dimension. We will show that although localization is seen in small world networks, to observe GPs, with 
power-law dynamics, we need networks of finite topological dimension. 



R2 


Figure 1 . Two lowest levels of the HMN2d hierarchical network construction. Dashed lines frame bottom level nodes, 
which are fully connected there, dotted lines frame the nodes of the next level. The solid lines denoted R1 are randomly 
chosen connections among the bottom level modules, ensuring single connectedness of the network, while those denoted R2 
provide random connections on the next level. Links can be directed. 

To generate a hierarchical modular network, we define l^ax levels on the same set of A = 4^'"“ nodes; on the level we 
define 4^ modules. We achieve this by splitting each module into four equal sized modules on the next level, as if they were 
embedded in a regular, two-dimensional (2d) lattice (see Figure 1). The probability pi that two nodes are connected by an 
edge follows pi ~ as in,^® where I is the greatest integer such that the two nodes are in the same module on level 1. After 
selecting the number of levels and nodes we fix the average node degree, b = (k) /2, and generate the adjacency matrix A by 
proceeding from the highest to the lowest level. We fill the submatrices with zeros and ones and copy them to appropriate 
diagonal locations of A. We allow unidirectional connections, which is more realistic. In fact, disorder associated with the 
link orientations turns out to be necessary to observe GPs. Finally, we connect the lowest-level modules with an edge, chosen 
randomly from nodes of li. This provides a short-linked base lattice that guarantees connectivity. 

In a preliminary study, the extra edges, corresponding to the base lattice, were not added. The resulting networks typically 
consist of a large number of isolated connected components; the GP effects observed in these structures are a consequence 
of fragmented domains of limited size. Measurements of axon-length distributions in several real neural networks show a 
large peak at short distances followed by a long flat tail. Thus there is a dense set of local edges in addition to a sparse 
network of long-range connections, best fit by an exponential function.^ Typical estimates indicate that ~ 90% of all cortical 
connections are formed at the local level (i.e., within a radius of 0.5mm); only the remainder leave the local gray matter and 
project to remote targets. Therefore, we also investigate a variant of HMN2d networks in which the lowest-level modules are 
fully connected, and there is a single link among the nodes of modules on level 12- To broaden further the range of structures, 
we study hierarchical modular trees (HMTs), which possess the minimum number of edges required for connectivity, and so 
have no loops. Construction of HMTs is described in the Supplementary material. 
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Relation to Benjamini-Berger networks 

Due to the embedding, there is a correspondence with Benjamini-Berger (BB) networks.^^ BB networks are generated on Eu¬ 
clidean lattices, with short links connecting nearest neighbors. In addition, the network contains long links, whose probability 
decays algebraically with Euclidean distance R: 

p{R)^R-r (2) 

Here we consider modified BB networks, in which the long links are added level by level, from top to bottom, as in.^® The 
levels : / = are numbered from the bottom to top. The size of domains, i.e., the number of nodes in a level, grows 

as Ni = 4^ in the 4-module construction, related to a tiling of the two dimensional base lattice. Due to the approximate 
distance-level relation, R ~ 2^ the long-link connection probability on level I is: 

E/=^(|) . (3) 

Here b is related to the average degree of a node (k). 

It is known that in a one-dimensional base lattice the BB construction results in a marginal case, i = 2, in which the 
topological dimension is finite and changes continuously with b. Eor i < 2 we have small world networks, while for s > 2 
the the topological dimension is the same as the base lattice (t/ = 1) in BB networks. The HMN-1 networks studied by 
Moretti and Munoz^^ are similar to the BB model, without the Id base lattice, but with the inclusion of a HMN topology. 
These authors simulated spreading models on HMN-1 with finite topological dimension at pi ^ and found GPs. Given 
the above mapping, this result is not very surprising, because this HMN corresponds to the s = 2 case, staying close to the 
percolation threshold of the network. To ensure connectivity, Moretti and Munoz added extra links to the large connected 
components;^'* they assert that the topological dimension remains finite, despite the additional links. In networks with finite 
sizes this is certainly true, however in the infinite size limit such a system is also infinite dimensional, so that we don’t expect 
true GPs. Eor the HMN-2, the number of connections between blocks at every level is a priori set to a constant value. In this 
case numerics of“^’^^ suggest GP effects. 

Here we embed the HMNs in 2d base lattices, which is closer to the real topology of cortical networks. In this case the 
threshold for marginality (i.e., continuously changing dimension and exponents) is expected to be at i = 4. We confirm this by 
explicit measurements of the topological dimension. Eurthermore, we study two-dimensional HMNs with different s values 
and find GP effects for finite topological dimensions, both at the percolation threshold (for s = 3), and for s — b, where the tail 
distribution of the link lengths decays very rapidly, in agreement with empirical results for axon lengths. 

In addition to the CP (m = 1), we also consider threshold models (m > 2), which are expected to be more realistic for 
neuronal systems. 

Dimension measurements 

To measure the dimension of the network we first compute the level structure by running a breadth-first search from several, 
randomly selected seeds. We count the number of nodes N{r) with chemical distance r or less from the seeds and calculate 
averages over the trials. We determined the dimension of the network from the scaling relation N ~ r'^, by fitting a power-law 
to the data for N{r). 

At s = 3 we observe a percolation threshold near {k) = 4, for which the topological dimension is finite: ~ 1.4 (see 
Eigure 2). Note that the curves with large (k) exhibit saturation, corresponding to a finite-size effect. The case s = 4 appears to 
correspond to a marginal dimension, for which, above the percolation threshold (k) > 6, one observes a continuously varying 
topological dimension (see Eigure 3). A more detailed, local slopes analysis via the effective topological dimension 


ln[A(r))/A(rO] 

ln(r/r') 


(4) 


is shown in the inset of Eigure 3, where r and r' are neighboring measurement points. Saturation to s-dependent dimensions 
can be read off the plot for r > 100. 

Random hierarchical trees also exhibit a finite topological dimension. Eigurere 4 shows N{r) for a set of ten independent 
random trees. The average of N{r) over the set of ten trees follows a power law to good approximation, with d ~ 0.72. (Note 
that in this case N{r) is an average over all nodes.) Eor regular trees, by contrast, N{r) grows exponentially with r, as shown 
in Eigure 5. Thus the regular tree construction results in a structure with infinite topological dimension. 
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Figure 2. Number of nodes within the chemical distance r in the s = 3 HMN2d models with l^ax = 9. Different curves 
correspond to different parameters: b= {k)/2as indicated. The dashed line shows a power-law ht for b = 2, with N{r) ~ r^-^, 
suggesting a topological dimension of t/ = 1.4 at the percolation threshold. 



Figure 3. Number of nodes within chemical distance r in HMN2d networks with s = 4 and / = 9 levels. Different curves 
correspond to different b = {k) /2 values as indicated. Inset: local slopes dgff of the N{r) curves, dehned in equation. 4. 


Dynamic simulations 

Contact process 

The contact process a Markov process dehned on a network, can be interpreted as a simple model of information 

spreading in social,epidemic spreading,^^^' or in brain networks.In the CP sites can be active or inactive. Active sites 
propagate the activity to their neighbors at rate X/k,or spontaneously deactivate with rate v = 1 (for simplicity). 

We perform extensive activity-decay simulations for HMN2d models with 2b — (k) = 4 and maximum levels: Imax = 
8,9,10. In these networks we use directed links between nodes, similar to real nervous systems. We follow p(f) in 10 — 100 
runs for each independent random network sample, starting from fully active initial states, and averaging p (f) over thousands 
of independent random network samples for each X. In the marginal long-link decay case at s = 4, a clear GP behavior is 
found (see Figure 6). 

We show that a GP can occur for more general parameters than those studied in,^‘* by following the density decay in 
networks with s = 3 and {k) = 4 (i.e., at the percolation threshold). Size-independent, nonuniversal power laws can be seen in 
Figure 7. 

When we increase the average degree (k), the GP shrinks to a smaller range of X values, as in,*^ tending toward a simple 
critical phase transition. However, it is hard to determine at precisely which value of (k) this happens. We note that the 
connectivity of the network is not assured, without the addition of extra links. Strictly speaking, however, with this extension 
i < 4 networks become inhnite-dimensional in the N limit. 

More importantly, GPs are also found in networks connected on the base level via short edges. We study the activity decay 
for s = 6, which corresponds to fast decaying tail distribution for the long links, preserving hnite topological dimension d. In 
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Figure 4. Number of nodes within chemical distance r in a set of ten random hierarchical trees with 262144 nodes (thin 
curves). The broad yellow curve is an average over the set of ten structures. 



r 

Figure 5. Number of nodes within chemical distance r in regular hierarchical trees with 1024, 4096,..., 262144 nodes 
(lower to upper curves). 


this construction the average degree at the bottom level is {k) ~ 6, decreasing systematically with I, so that (k) ~ 1 for Imax- 
Note that that the ratio of short to long links is ^ 0.11, in agreement with results for real neural networks.^^ Simulations again 
yield size-independent power-law decay curves, confirming GP behavior, as shown in Figure 8. 

We also determined the effective decay exponent in the usual manner (see^^), via 


ln[p(f)/p(f')] 




ln(f/f') 


(5) 


where p(f) and p(f') are neighboring data points. The critical point can be located at Xc = 2.53(1), showing mean-held scaling: 
p (f) ^ f^*. Above this threshold power-laws can still be seen for smaller sizes Umax = 8,9), up to A ~ 2.55, but corrections 
to scaling become stronger; the curves for Imax = 10 exhibit saturation. This is surprising, because in a disordered system 
with short ranged interactions one would expect a critical phase transition with an ultra-slow, logarithmic scaling. However, 
recent studies of the CP in higher dimensions hnd mean-held criticality and GP.^^ Our result suggests that the topological 
heterogeneity generated by the long edges is not strong enough to induce an inhnite-randomness hxed point. Otherwise we 
must assume very strong corrections to scaling in this case. 

Disorder due to the randomly chosen orientations of the links turns out to be relevant. For symmetric links our simulations 
yield nonuniversal power laws, which appear to be sensitive to the system size, suggesting the lack of a true GP in the inhnite- 
size limit. 


Burstyness in the CP 

We study the distribution of inter-event times in the CP on asymmetric, HMN2d networks with s = 6, in a manner similar to 
that described in.^'^ Starting from fully active states, on networks of size Imax = 10, we measured (At) between successive 
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Figure 6. CP on asymmetric HMN2d networks with i = 4; decay of activity for X values as indicated. System sizes 
Imax = 8,9,10 (thin, medium, and thick lines, respectively). Size-independent power laws are evidence of a Griffiths phase. 



Figure 7. CP on asymmetric HMN2d networks as in Figure 6, but for s = 3. 

activation attempts. We followed the evolution up to tmax = 2^^ Monte Carlo steps (MCs), averaging over 1000 — 2000 
independent random networks with 10 — 100 runs for each. Througout this paper time is measured by MCs. As Figure 9 
shows, fat-tailed distributions, P{At) ^ emerge in the GP for 2.46 < A < 2.6, while P{At) decays exponentially 

outside this region. The slopes of the decay curves are determined via least-squares fits in the window 20 < Af < 7000. For 
X —2.5 we find x = 1.753(4), while for A =2.52 the exponent is slightly larger; x = 1.81 (1). These values are close to the 
auto-correlation exponent of the critical 1-1-1 dimensional CP as in,^'* but exhibit deviations due to the heterogeneities. This 
is understandable, since the effective dimension of this HMN2d is close to one. The scaling variable P*{Al) = (Af)*-^®P(Af) 
exhibits log-periodic oscillations. This is the consequence of the modular structure of the network. Furthermore, as in other 
GP models, logarithmic corrections to scaling are expected. 

We expect similar nonuniversal, control parameter dependent tails in P{At) in the GPs exhibited by the other HMN2d 
networks. Furthermore, as shown in,^'* power-law distributions should arise for other initial conditions, such as localized 
activity. This suggests that bursty inter-communication events in brain dynamics arise spontaneously near the critical point, in 
the GP. 

CP on random hierarchical trees 

We simulate the CP with symmetric links on random hierarchical trees (RHTs) of 262 144 nodes. We first perform quasista¬ 
tionary (QS) simulations,^^ of the CP on a single RHT. For one structure this yields Ac — 2.72; for another, independently 
generated structure of the same kind, we find Ac ~ 2.76. 

Our principal interest is in the decay of activity starting from all nodes active. The decay of activity on a single RHT 
appears to follow the scenario familiar from the CP on regular lattices: power-law decay is observed at Ac but not at nearby 
values, as illustrated in Figure 10. The randomness associated with a single RHT appears to be insufficient to generate a 
Griffiths phase. 
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Figure 8. CP on asymmetric HMN2d networks with s = 6 and (k) ~ 4; decay of activity for X values as indicated. System 
sizes as in Figure 6. Size-independent power laws are observed for 2.45 < A < 2 .53. Inset; local slopes of the curves in the 
main plot. In the GP Ueff tends to nonuniversal values with logarithmic corTections. At Xc = 2.53(1) we observe convergence 
to Ueff = 1; above this threshold the effective exponents of larger systems veer upwards toward zero as t tends to infinity. 



Figure 9. CP on asymmetric HMN2d networks with s = 6. Main plot: probability distribution, P{At), of inter-event times 
for X values as indicated; system size l^ax = 10. Power-law tails are evident for 2.46 < A < 2.6, with continuously changing 
exponents. Right inset; local slopes of the same curves, defined similarly as in equation (5) for; A = 2.46,2.48,2.5,2.52,2.7. 
Left inset; P*{At) = (Af)' ™P(Af) for A = 2.5,2.51,2.52. 


By contrast, evidence of a GP is found if we average over many RHTs. The activity averaged over a large set (10^ - 10^) 
of independent realizations, each on a different RHT, shown in Figure 11, decays asymptotically as a power-law over a fairly 
broad range of subcritical A values. The decay exponent a extrapolates to zero at A = 2.760(2) ~ Ac, as shown in the inset. 
We note that the average is over all realizations, including those that become inactive before the maximum time of 2 x 10® 
MCs. If we instead restrict the average of p(f) to trials that survive to time f (or greater), the result is a stretched exponential, 
p (f) ^ exp[—Cf^], where C is a constant and the exponent varies with A. For A = 2.70, for example, ~ 0.25. 

Threshold model simulations 

Threshold models represent an attempt to capture the activation threshold of real neurons, by requiring at least m active 
neighbors associated with incoming links to activate a node with probability A. In case of active nodes spontanous deactivation 
occurs with probability v written in the reaction-diffusion model notation as; 

(a) niA —>■ (m+l)A with probability A, 

(b) A —> 0 with probability v. 

We use SCA updating, as in.^® Since there is no spatial anisotropy (which might generate activity currents), we can assume that 
the SCA follow the same asymptotic dynamics as the corresponding model with random sequential updating. Thanks to the 
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Figure 10. CP on a random hierarchical tree; activity versus time for X — 2.73, 2.74,..., 2.77 (lower to upper); system size: 
N = 262144. Power-law decay is evident only for A = 2.75, close to the estimate A^ = 2.76 derived from QS simulations. 
Each curve represents an average over 100 realizations; all realizations are performed on the same network. 


synchronous updates, we can speed up the simulations by ~ Up times by distributing the nodes among np multiprocessor cores; 
by swapping the random number generation on parallel running graphics cards we obtain a further reduction of 50% decrease 
in the simulation times. The spatio-temporal evolution of the HMN2d network with m = 6, in a single network realization, 
starting from a fully active state, is shown in Figure 12. After a sharp initial drop in activity, due to spontaneous deactivation 
of nodes with few neighbors, one observes domains (modules) on which activity survives for a long time, suggesting rare 
region effects. 

Results for threshold models 

We begin by discussing the m = 2 case, for which extensive simulations are performed. We generate networks with average 
degree {k) = 24 and s = 6. Note that in this case values of (k) higher than the percolation threshold are needed to avoid 
modules having separated activities. As before, we averaged over hundreds of independent random networks and thousands 
of independent realizations. We followed the density up to t^ax = 10^ MCs, starting from a fully active initial condition. 

In this case the mean activity density decays more rapidly than in mean-field theory, and is size-independent (see Figure 1 3). 
For large branching probabilities, (p(f)) seems to take a constant value, suggesting an active steady state, but at late times 
some deviation is observed, possibly due to finite-size effects. 

Homogeneous triplet creation models (i.e., m = 3) are expected to exhibit a mean-field-like discontinuous phase transition 
in two or more dimensions(see for example^^). Disorder induces rounding effects, producing continuous phase transitions 
or GPs.^^’^^ We study a threshold model with m = 3 and s = 6, using (k) = 36. Our results are similar to those for m = 2: 
non-universal power-laws, again suggesting a GP. 

Quenched Mean-field approximation for SIS 

Heterogeneous mean-field theory provides a good description of network models when fluctuations are irrelevant. This ap¬ 
proximation is attractive because is can be solved analytically in many cases; the results agree with simulation.This 
analysis treats nodes with different degrees as distinct, but finally averages over all degree values, providing a homogeneous 
solution for the order parameter. To describe the effects of quasi-static heterogeneities in a more precise way, the so-called 
Quenched Mean-Field (QMF) approximation was introduced.For SIS models this leads to a spectral analysis of the adja¬ 
cency matrix j of the network. The susceptible-infected-susceptible (SIS)®** model is similar to the CP, except that branching 
rates are not normalized by k, leading to symmetric governing equations. A rate equation for the SIS model can be set up for 
the vector of activity probabilities p, (f) of node i at time f; 

= -Pi{t) + A(1 - P;(f)) E ■ (6) 

Here the Wij = wji are weights attributed to the edges. For large times the SIS model evolves to a steady state, with an order 
parameter p = (p,). Since this equation is symmetric under the exchange i j, a spectral decomposition can be performed 
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Figure 11 . CP on a random hierarchical trees; activity versus time for A = 2.60, 2.62, 2.64, 2.68, 2.69,..., 2.75 (lower to 
upper); system size N = 262144. Each curve represents an average over 10^ - 10^ realizations, each performed on a different 
network. Inset; decay exponent — a versus A. 


on a basis of orthonormal eigenvectors e{y). 


Pi = ^c(A)e;(A). 

A 

Non-negativity of the matrix Bij = AtjWij assures a unique, real, non-negative largest eigenvalue yM- 
In the long-time limit the system evolves into a steady state and we can express the solution via Bjj as 


(7) 




^ L; BjjPj 
1 + 


( 8 ) 


Stability analysis shows that p, > 0 above a threshold Ac, with an order parameter p = {pi). In the eigenvector basis equation (8) 
can be expanded by the coefficients c(A) as 


:(A) = A EA'c(A') f- ^ 


(9) 


and near the threshold we can express p, ^ c,(Ai )e,(Ai) with the principal eigenvector. In the QMF approximation 1 /Ac = y\ 
and the order parameter can be approximated by the eigenvectors of the largest eigenvalues 


p(A) ~ ci\A-\-(22A -f... , 
where A = AAi —1^1 with the coefficients 


( 10 ) 


( 11 ) 


(=1 


(=1 


As proposed in,^® and tested on several SIS network models, localization in the active steady state can be quantified by 
calculating the inverse participation ratio (IPR) of the principal eigenvector, related to the eigenvector of the largest eigenvalue 
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Figure 12. Spatio-temporal evolution of the density of the m — 6 threshold model, proportional to the color coding on the 
right. The simulation is started from an active state, with v = 0.9 and X = 1. 



Figure 13. Decay of activity in the m = 2 threshold model with s = 6, and (k) = 24. The curves correspond to branching 
rates: X = 0.65,0.66,0.67,0.68,0.83 (lower to upper) and v = 0.7 fixed. Levels: Imax = 8,9 (thin, thick lines). 
Size-independent power-laws, reflecting a GP are observed. 


e(yi) of the adjacency matrix as 

( 12 ) 

i=l 

This quantity vanishes ^ 1 /N in the case of homogeneous eigenvector components, but remains finite as oo if activity is 
concentrated on a hnite number of nodes. Numerical evidence has also been provided for the relation of localization to RR 
effects, slowing down the dynamics to follow power-laws}'' 

We study localization of the SIS in the HMN2d models introduced in the previous section*s. We analyze the eigenvectors 
of Bij for b — 2{{k) ~ 4), varying s, using system sizes ranging from N = 256 to N — 262144. In particular, we performed a 
FSS scaling study of the IPR in these systems. First we determined the behavior on the small world network of^^ corresponding 
to s = 3/2. As one can see in Figure 14, the localization at small sizes disappears as I{N) ~ 1 /A^. By contrast, for s = 3 
and s = 4, the graphs are hnite-dimensional, and the IPR increases with N, tending to a hnite limiting value, suggesting 
localization. 

One might question the relevance of network models with small connectivity to mammalian brains, in which {k) is on the 
order of 10^. To answer this we study s < 4 models with higher connectivity, i.e., (k) ~ 50. As Figure 15 shows, the sign of 
localization, which is weak but nonzero for (k) = 4, now disappears. Next, we add random weights Wij, distributed uniformly 
over the interval (0,1), to the links of the networks. A consequence of this explicit disorder is localization even in highly 
connected networks, as shown in Figure 15. This result is in agreement with the expectation of limited sustained activity in 
brain networks,^® meaning that the link disorder prevents over-excitation of a network of high connectivity. 
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Figure 14. Finite size scaling of the inverse participation ratio in weakly coupled HMN2d models, with maximum levels 
Imax = 4,5, ..9. The s = 3 (bullets) and i = 4 (boxes) results suggest localization (finite I) in the infinite-size limit. Lines are 
power-law fits to the data. For i = 1.5, corresponding to the symmetrized, small-world network (model-6 of^^) no evidence 
of localization is seen. 



Figure 15. Finite size scaling of the inverse participation ratio in HMN2d models with higher average degree {{k) ~ 50) for 
maximum levels: Imax = 4,5,6,7,8. Bullets: s = 3 with uniform randomly distributed weights; boxes: without weights. 
Diamonds: i = 4 with randomly distributed weights. Lines show power-law fits to the data. In the unweighted case, no 
localization effect can be seen, and I decays linearly with N. 


Localization suggests rare-region effects, thus a dynamic GP. Nevertheless, simulations of the CP on such networks do 
not show extended regions of power-laws for s < 4, but rather a nontrivial (non-mean-field) continuous phase transition 
(Figure 16). Decay simulations for I = 9 yield p(t) ~ f-0-50{i) 2.588(1), albeit with a cutoff due to the finite system 

size. This agrees with our result for ID BB networks with s = 3.^^ Spreading simulations starting from single, randomly 
placed seed result in the survival probability scaling at this critical point: P(f) ~ f-0 50(1) symmetry a = 5 holds 

to within uncertainty). Here the density increases initially as p(f) ~ f0.23(i) Qjjg understand these results, noting that 
the localization values are rather small, I ^ 0.08, so that RR effects are too weak to generate observable GP effects in the 
dynamics. 

Naturally, for i > 4 networks, where already the topological disorder is strong enough to create a GP, inhomogeneities in 
the interaction strengths amplify the RR effects. This induces GPs even for larger values of (k), as well as in the threshold 
models. Further studies of disorder effects are under way. 

Conclusions 

We investigate the dynamical behavior of spreading models on hierarchical modular networks embedded in two-dimensional 
space, with long links whose probability decays as a power-law with distance. This corresponds to an exponentially decreasing 
connection probability, as a function of the level in the hierarchical construction. The aim of this study is to understand the 
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Figure 16. Evolution of the activity in the weighted CP on HMN2d networks with s = 4,1 — 9, and k = (50), for 
X = 2.585, 2.588, and 2.59 (lower to upper). Main plot; decay in case of a fully active initial state. Inset; growth of activity 
starting from a single active node for the same values of X. Dashed lines are power-law fits for Xc — 2.588. 


effects of intrinsic and topological disorder, in particular, extended critical regions in the parameter space, without any (self) 
tuning mechanism. 

If we eliminate the underlying lattice structure we observe power-law dynamics for networks near the percolation threshold, 
when the effective dimension is hnite. However, size-independent power-laws, corresponding to GPs, are only seen if we have 
directed links. Since connectivity of these structures is not guaranteed, we also study random hierarchical trees, with full 
connectivity. GPs are observed upon averaging over many independent network realizations. The relation to brain networks 
can be envisaged in the large-size limit, if we regard independent realizations as (almost decoupled) sub-modules of the entire 
brain. 

When we ensure connectivity via short edges on the lowest level, we find GPs for rapidly decaying long-link probabilities. 
Both of these network assumptions are in accordance with empirical brain networks. Above the GP, at the critical point, we 
find mean-field-like decay of activity, in agreement with results on the CP with quenched disorder on higher-dimensional, 
regular lattices. We have also shown that bursty behavior arises naturally in the GPs, due to autocorrelations which decay via 
a power-law. 

We perform a quenched mean-held study of the SIS model on these networks; in the SIS, nodes are connected symmet¬ 
rically. Finite size scaling of the inverse participation ratio suggests localization on large-world networks and de-localization 
on small world structures. However, when we add intrinsic weight disorder, localization can be seen even on small-world 
networks. Weight disorder is again to be expected in real brain networks, since the strength of couplings varies over many 
orders of magnitude. Despite this, we saw no GP effects in the dynamics of weighted CPs with i = 4. Instead, we hnd a 
nontrivial critical scaling, as has already been observed in other networks. The possibility of a narrow GP in this case is an 
open issue. 

In conclusion, we believe our synthetic HMN2ds are closer to experimental brain networks than previously proposed 
models, and hnd numerical evidence for GPs in extended phases in simple models with spreading dynamics. Although we 
eliminate any self-tuning mechanisms, we still hnd nontrivial slow dynamics as well as localization of activity, which is crucial 
for understanding real brain network data. 
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Supplementary information 

Although the set of nodes is isomorphic to a square lattice of N = 2^x2^ sites, we shall label the nodes not by their Cartesian 
coordinates but rather using a base-4 notation. At level I, the full set of 4^ nodes is divided into four quadrants labeled 0, 
1, 2, and 3, proceeding counterclockwise from the lower left. At level / — 1, each of the quadrants is similarly divided into 
four subquadrants, labeled in the same manner. Division into subquadrants continues all the way down to level 1, in which 
the elements are individual nodes. A given site can be specified by the labels of the quadrant, subquadrant, etc. to which it 
belongs. Denoting the labels at levels 1,2,...,L, by gi, ...,gL, the address of a site is n = g/4^^* (see Lig. S17). 


15 

14 

11 

10 

12 

13 

8 

9 

3 

2 

7 

6 

0 

1 

4 

5 


Supplementary Figure S 17. Definition of node labels for a network with L — 2 

Consider a set of four nodes, 0, 1, 2, and 3. There are six possible links between them; (0,1), (0,2), (0,3), (1,2), (1,3), and 
(2,3). A tree with four nodes has exactly three links. Of the ( 3 ) = 20 subsets of three links, sixteen yield a connected tree 
[for example: {(0,1), (1,2), (2,3)}], while four correspond to a cycle of three nodes, leaving the fourth node isolated [example: 
{(1,2), (1,3), (2,3)}]. To construct a tree of four nodes, we choose one of the sixteen link sets at random. 

A hierarchical tree is constructed by first linking the four quadrants via three edges. The link set iMj is chosen at random. 
Then, for each link, we choose sites at random within each of the two quadrants connected by the link, to serve as the connected 
nodes. Now we repeat this process within each of the subquadrants, and so on, until we reach the basic modules of four sites. 
The latter are again connected by sets of three links, chosen independently from the basic collection of sixteen link sets. At 
the end, each basic module is connected internally, and to a module at the next level, etc., so that we have a connected graph 
of N nodes and N —1 edges. Although we have chosen to begin the construction at level L, we could equally have begun at 
level 1, or some intermediate level: the choices of link sets and nodes at the various levels are mutually independent. These 
steps are illustrated in the Ligs. S18, S19, S20. 
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Supplementary Figure S 18. First step in constructing a network with L = 2: define the link set at level 2. 



Supplementary Figure S 19. Second step: choose nodes consistent with link set at level 2. At the next step (not shown) 
we choose link sets for each of the four-site modules. 



step: choose nodes at level 1 consistent with the corresponding link sets. 


At level 1 the number of links per node is 3/4; at level j, there are 4^^-' blocks connected via — 1 links. The number 
of links per node at this level is therefore (4^^^ — l)/4^ ~ 1/4T Let pj denote the probability that a randomly chosen edge 
linking blocks at level j be present. At level 1, links connect nodes within four-site modules. Since there are six possible links, 
and since just three are present within each module, we have pi — I /2. At level 2, a link connects nodes within a 16-node 
module; at this level a node is linked to one of the 12 nodes outside its basic 4-site module. There are (16 • 12)/2 = 96 possible 
links, and since only three are present, P 2 = ^ /32. At level j, an edge links nodes within a module containing 4-' nodes. The 
number of possible links at this level is 4'' (3 ■4''^*)/2, so that pj = 1/2"'^-'^^, and, in general, pj/pj-i = 1/16. Since a tree 
represents a minimally connected structure, pj cannot decay faster than this rate, in any connected hierarchical network based 
on four-node modules. 
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Although our principal interest is in random trees, as described above, we may also consider regular trees; these are 
constructed by using the same link set at each level. The resulting structure turns out to have inhnite topological dimension, 
as discussed in Section ’’Dimension measurements”. 
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